Spatial analysis in R

Workflow

After several tentative of connecting Google Earth Engine and R via rgee, I had to abandon this idea and make a different workflow. The initial idea was to export geojson files from GEE and plot them in R. In some cases there were errors with GEOMETRYCOLLECTION features that blocked both leaflet and mapview. I then switched to this workflow, using QGIS for intermediate calculations.

flowchart TB
  A[Dataset] --> B(GEE)
  B --> C[export tif]
  C --> D[QGIS]
  D --> E[Zonal statistics]
  E --> F(import in R)
  F --> G[plot with leaflet]

Printing vector data

The first example will use a geojson file created from Google Earth Engine indicating the population in each administrative district of Sweden.

Printing with ‘sf’ package

plot(data_map[,14],max.plot = 14)

Printing with ‘ggplot’ package

ggplot(data_map) + geom_sf(aes(fill = sum)) +
  scale_fill_viridis() + theme_bw()

Printing with ‘leaflet’ package

leaflet()%>%
  addTiles()%>%
  addPolygons(data= zonal)

Printing with ‘mapview’ package

#mapview(zonal, zcol = "_mean")
#mapview(map, zcol = "count")
#mapview(pop, zcol = "sum")
#mapview(dens, zcol = "mean")
#mapview(dens_multi, zcol = "mean")

Both leaflet and mapview struggle with GEOMETRYCOLLECTION features.

Printing raster

r <- rast("soilT_swe_1000.tif")
plot(r)

Plot using terra

map <- rnaturalearth::ne_states("Sweden", returnclass = "sf")
sextent <- terra::ext(map)
r <- terra::crop(r, sextent)
r <- terra::mask(r, vect(map))
plot(r)

plot(map)

Extracting raster values at point

v <- vect("admin.shp")
points <- crds(centroids(v))
plot(r)
plot(v, add = TRUE)
points(points)

# data frame with the coordinates
points <- as.data.frame(points)
valuesatpoints <- extract(r, points)

Averaging values

df_area <- extract(r, v, na.rm = TRUE)

# Average raster values by polygon
v$avg <- extract(r, v, mean, na.rm = TRUE)$SoilTemp10_40cm_tavg_count

# Area-weighted average raster values by polygon (weights = TRUE)
v$weightedavg <- extract(r, v, mean, na.rm = TRUE,
                 weights = TRUE)$SoilTemp10_40cm_tavg_count

library(ggplot2)
library(tidyterra)

# Plot average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = avg)) +
  scale_fill_terrain_c()

# Plot area-weighted average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = weightedavg)) +
  scale_fill_terrain_c()

Printing raster with leaflet

rb <- raster::brick(r)

pal <- colorNumeric("YlOrRd", values(r),
                    na.color = "transparent")
leaflet() %>% addTiles() %>%
  addLegend(pal = pal, values = values(r), title = "elevation")%>%
  addPolygons(data = zonal$geometry,
              weight=1,
              group = "Admin")%>%
    addRasterImage(rb, colors = pal, opacity = 0.4, group="Temp") %>%
   addLayersControl(
    overlayGroups = c("Temp", "Admin"),
    options = layersControlOptions(collapsed = FALSE)
  )

https://rstudio.github.io/leaflet/articles/showhide.html